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ABSTRACT 

We consider the local dynamics of a realistic neutron star core, including composition 
gradients, superfluidity and thermal effects. The main focus is on the gravity g-modes, 
which are supported by composition stratihcation and thermal gradients. We derive 
the equations that govern this problem in full detail, paying particular attention to 
the input that needs to be provided through the equation of state and distinguishing 
between normal and superfluid regions. The analysis highlights a number of key issues 
that should be kept in mind whenever equation of state data is compiled from nuclear 
physics for use in neutron star calculations. We provide explicit results for a partic¬ 
ular stellar model and a specific nucleonic equation of state, making use of cooling 
simulations to show how the local wave spectrum evolves as the star ages. Our results 
show that the composition gradient is effectively dominated by the muons whenever 
they are present. When the star cools below the superfluid transition, the support for 
g-modes at lower densities (where there are no muons) is entirely thermal. We confirm 
the recent suggestion that the g-modes in this region may be unstable, but our results 
indicate that this instability will be weak and would only be present for a brief period 
of the star’s life. Our analysis accounts for the presence of thermal excitations encoded 
in entrainment between the entropy and the superfluid component. Finally, we discuss 
the complete spectrum, including the normal sound waves and, in superfluid regions, 
the second sound. 
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1 INTRODUCTION 


It is well established from solar physics and the hugely successful Soho mission that helioseismology ( Garcia et al.|20 T^ may 
provide insights into the physics beneath a star’s surface. Following on from this, missions like Corot and Kepler | Aerts|2015[ ) 
have firmly established asteroseismology as a precision science. However, in order to make use of the seismology strategy one 
must have a good theoretical understanding of the physics involved and the nature of the star’s various oscillation modes. One 
also relies on nature to provide a mechanism that excites these modes to a detectable amplitude. In the case of the Sun, the 
interior physics is (by now) relatively well understood and we know that interior dynamics can be excited by convection. The 
situation is quite different when it comes to extreme objects like neutron stars. These, highly degenerate, systems also support 
a plethora of oscillation modes and it has been demonstrated how different aspects of supranuclear physics may influence 
these stars’ oscillations ( McDermott et ^|1988 Kruger et al.|[2M4 l. However, we are still quite far from a complete picture 
(partly because the physics of the deep core of these stars remains poorly understood and partly because of difficulties in 
building models that account for the rich physics we know we have to include). 

The neutron-star problem is interesting (and topical) since an oscillating neutron star may radiate gravitational waves 
at a level that would be within reach of the advanced ground based detectors that are now coming online. In practice, this 
would require large-amplitude oscillations due to some kind of instability being active ( Andersson et al.|2011 |. Given this, it 
is understandable that much recent effort has been focussed on various neutron-star instabilities (e.g. the gravitational-wave 
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driven instability of the inertial r-modes (Andersson & Kokkotas 2001 Ho et al. [2011 Alford & Schwenzer 20141 or the 
f-modes ( Passamonti et al.|2013 Doneva et al. 12013 1). 

The matter composition in the star plays a key role in determining the fluid dynamics. Basically, the richer the “chemistry”, 
the more complex the phenomenology may be. Variations in composition generally enables moving fluid elements to experience 
buoyancy whenever they are out of equilibrium with their surroundings. This leads to the presence of the so-called gravity 
g-modes. In more familiar settings, e.g. the Sun, the main stratification is thermal. In contrast, neutron stars tend to be so cold 
that thermal effects are irrelevant. Nevertheless, local fluid elements may experience buoyancy due to composition gradients 
| |Reiseneggerj^^^Goldreich| 19^ I. This problem has been studied at some level of detail for simple neutron star models (see, for 
example, Passamontret*^.|2009 Gaertig Kokkotas|2009 1 and the relevance for various astrophysical scenarios, ranging from 
hot proto-neutron stars ( Ferrari^^et^^alj20^ Burgioet al.|20lT \ and adolescent stars ( Kruger^et^^lj20^ I, to tidal interactions 
( Kokkotas fc Schafer|1995 Flanagan fc Racine|2007' I, various instabilitie s (|Lai|1999| |Weinberg et al.|2013[ ), and the dynamics 
associated with core collapse and the generation of gravitational waves | Ott et al.|2006 1, has been considered. More realistic 
neutron-star models have not yet been considered in particular detail, although it is known that interfaces associated with 
phase transitions will lead to the presence of a family of modes closely related to the g-modes ( Miniutti et al.|[2003 l. It has 
also been demonstrated that the onset of superfluidity has a key influence on the buoyancy that supports the g-modes. In 
simple models, the support for the g-modes may, in fact, disappear altogether | Lee|199^ Andersson fc Comer|200T |. 

The traditional view is that one would not expect the convection that drives observed oscillations in main-sequence stars 
like the Sun to operate in a mature neutron star. A neutron star becomes stably stratified almost immediately (after 100s or 
so) after birth and, even though the composition will affect the dynamics, the various oscillation modes are expected to be 
convectively stable. This conclusion appears to be challenged by recent work (Gusakov fc Kantor||2013 Kantor & Gusakov 


|2014p that suggests that the presence of superfluidity in the neutron star core may lead to local regions becoming convectively 
unstable. At first sight, this result is surprising, so it seems important to establish (first of all) whether it is correct. If the 
presence of superfluidity can, indeed, trigger a convection phase then we will need to understand the implications this may 
have for (say) the star’s thermal evolution and other observable phenomena. The present investigation aims to address the 
hrst of these questions. 

We consider the local dynamics of a realistic (outer) neutron star core, focussing on the buoyancy experienced by fluid 
elements and the associated gravity g-modes. Our model accounts for the presence of neutrons, protons, electrons and muons, 
and we include finite temperature effects by introducing an entropy component. However, we assume that the dynamics we 
are considering are sufficiently fast that thermal conductivity can be ignored. In essence, this means that i) we can decouple 
the slow evolution associated with cooling, and ii) if we allow the neutrons to becomes superfiuid then the entropy is locked 
to any normal fluid in the mixture. Since the protons are expected to form a superconductor (we neglect all electromagnetic 
aspects here) in the regions we are considering, the entropy will be carried by the electrons/muons. Once we have developed 
the detailed framework for analysing the local (plane-wave) dynamics of this model, we provide detailed results that combine 
an actual cooling simulation with our derived dispersion relation. This should allow us to establish whether it is realistic to 
expect that convectively unstable regions may, indeed, exist at some point of a neutron star’s early life. 

The motivation for this work is to understand i) the effect that the interior composition has on the local dynamics of 
the fluid in the star, and ii) how the onset of superfluidity impacts on the results. However, there are a number of crucial 
technical aspects to the analysis. Our discussion highlights the microphysics information that is required to study this problem 
in particular, and model detailed neutron star dynamics in general. The implementation of this information may be somewhat 
technical (and hence it is mainly discussed in Appendices), but it is important to understand to what extent available 
neutron-star equation of state models provide the required information. 

Throughout the discussion, we work in a coordinate basis, expressing vector quantities in terms of their components. 
We denote space-time indices by lowercase italics starting from the beginning of the alphabet, a, &, c,.... Spatial indices are 
also lowercase italics, but start from In both cases the Einstein summation convention is assumed. The signature 

of the metric gab is ,-I-,-I-,-I-], and the covariant derivative associated with this metric is denoted by Va. Different fluid 
components are distinguished by means of a constituent index; a Roman letter x, y,.... Specifically, we use n for neutrons, p 
of protons, e for electrons, g for muons and s for entropy. These matter indices are not summed over when repeated. 


2 KEY INGREDIENTS OF THE MODEL 
2.1 Relativistic multifluid dynamics 

As we want to allow for the presence of superfiuid components, it makes sense to build our model within the general framework 
for relativistic multi-fluid systems (for a review see | Andersson fc Gomer|2007p . This means that we take as our starting point 
the convective variational approach to relativistic fluids in which the key variables are the various fluxes n“, where x is a label 
that identihes the fluid (in the following it will be n, p, e, /i or s). In general, these fluxes do not have to be conserved, and 
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we would have 

Va< = 7x , (1) 

where 7x follows from the relevant reaction rates. The problem simplifies in two extreme limits. In the first limit, the reactions 
are much faster than the dynamics. This means that a moving fluid element has time to equilibrate to its surroundings 
and hence its composition changes accordingly during the motion. In the opposite limit, when reactions are slow, the fluid 
element retains its chemical identity. In this latter case, the mismatch between the composition of the fluid element and the 
neighbourhood leads to buoyancy. This can either serve to push the fluid element back in the direction it came from, which 
leads to a stable oscillation, or push it on in the direction of travel, which leads to an instability and large-scale convection. 
We will assume that we are working in the slow-reaction limit, and so ignore the 7x and take the individual fluxes to be 
conserved. 

In the text-book approach to relativistic fluid dynamics one assumes that friction locks the different species in the fluid 
together, and as a result they move with a common four-velocity However, iu a maturing neutron star the neutrons will 

become superfluid. This means that they will not experience friction (at least not in the usual sense) and hence they can 
drift relative to the other components. In principle, the same is true for the protons, which will become superconducting. 
However, in order to simplify the problem one may assume that the charged components are locked electromagnetically, 
forming a charge-neutral conglomerate. In essence, this means that we ignore the electromagnetic field (as there will be no 
charge currents). This assumption may not be fully justified (as neutron stars have magnetic fields and so must support charge 
currents), but we will nevertheless make it as it keeps the problem tractable. 

Assuming that protons, electron, muons and entropy co-move while neutrons are superfiuid, we have 

= nxu“ , where x = p, e, /r, s , (2) 

and 


= (m“ -b u“) 


(3) 


where is the relative velocity of the neutrons {u°'Va = 0). Here and in the following we assume that this relative velocity 
is small enough that we can ignore red-shift factors in these expressions (u^ <C 1). This should be true for all situations 
of astrophysical interest. The number densities Ux are then (effectively) all determined by the same observer. Above the 
superfluid transition, the neutrons are obviously locked to the other components. 

The momenta that are conjugate to the particle fluxes follow from an energy S (which in turn is determined by the 
equation of state). This leads to 




X 

a 



y 7^ X . 


(4) 


In general, this allows for the presence of entrainment; an effect that describes how one fluid is dragged along as another 
fluid moves. In a neutron star core, the protons and neutrons will be entrained due to the strong interaction. Later, we will 
demonstrate that the multi-fluid formalism accounts for this effect in a very intuitive way. However, in order to keep the initial 
analysis tractable, we will not account for the entrainment at this point. This means that we have the momenta 


= MxMa , X = p,e,/r, s , 


(5) 


where /ix are the respective chemical potentials and it is worth noting that the temperature is T = /is. We also have, when 
the neutrons are superfluid, 

fJ-a = /in {Ua + Va) ■ (6) 


In the absence of mechanisms coupling the components, the equations of motion for this system are represented by the 
conservation laws 0 (although with 7 x = 0) and the individual momentum equations; 

/“ = 2<V[,/i"6] =0 , (7) 

where the square brackets denote anti-symmetrisation. However, when some mechanism couples various components one has 
to account for the relevant force (/x yf 0). Alternatively, one can appeal to Newton’s third law and work with relevant 
combinations of the equations (such that the equal and opposite coupling forces cancel). In the problem we consider here, 
the upshot of this is that it is natural to work with the total momentum equation (obtained by adding all the component 
equations) and the momentum equatiou for the superfiuid neutrous (where we ignore coupling mechanisms like the vortex 
mediated mutual friction ( Mendell|1991 Andersson et al.]|2006 l). The first of these equations contains the same information 
as the divergence of the stress-energy tensor. 

Above the superfiuid transition temperature one might as well work with 
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where 

T“'’ = (p + e)uV+p^7“^ (9) 

where p is the pressure and e is the energy density. In equilibrium, the equation of state for matter is given by a barotrope 
p = p(e), but as soon as we account for perturbations we need to account for the detailed composition. This is an important 
point, which we will return to later. 

Below the superfluid transition, we can still use the sum of the momentum equations (although it will take a slightly 
different form as there will now be explicit terms of order u“). In order to describe the new degree of freedom associated with 
the superfluid neutrons we complement this system by /“ =0. 

We also have a set of thermodynamical identities. As we are neglecting terms of second order in the energy density 
and pressure obey the usual Gibbs relation (the integrated first law of thermodynamics) 

£ + P = U-x/Tx . (10) 

X 

From thermodynamics principles (cf. the definition of the conjugate momenta and the chemical potentials) we have 

d£ = ^Mxdnx, (11) 

X 

which means that 

dp = ^nxd/rx. (12) 

X 

These relations lay the foundation for our analysis. 


2.2 Equilibrium configurations 


Let us now turn to the problem of building a relativistic star. We want to consider the perturbations of a non-rotating star 
in dynamical and thermodynamical equilibrium (at least on the timescale of the dynamics we are considering). That is, our 
background model is a spherical star in which all fluids co-move (w“ = 0). Considering the general metric for a spherical star 
(taken to be fixed later, meaning that we work in the Cowling approximatiorQ ; 

ds^ = —-I- + P [dO^ + sin^ Odp^'j , (13) 

it is easy to see that the background 4-velocity and the conjugate momenta are: 

pi = ■ (14) 

Meanwhile, the total pressure is obtained from 

p'= - (p + e)^'= - (p + e)p- (15) 


Here, and in the following, primes denotes radial derivatives. We have also introduced the gravitational acceleration, g, for 
later convenience. Combined with the definition of A in terms of the mass enclosed within radius r; 

Mr) = ^ (l - ’ (16) 

and 

^ (m -I- 4npP) , (17) 

these are the familiar Tolman-Oppenheimer-Volkoff equations. For later convenience, it is also worth noting that the Einstein 
equations lead to 


Finally, if the neutrons are superfluid, 

/ 

Mn 


A' -b <&' = 


(p + e) • 


then we also have ( Gusakov fc Andersson|2006 ) 

= —pn^' - Pn = pnP = COUStant . 


(18) 

(19) 


^ It is worth noting that we should possibly be a little bit careful here. As shown by |Finn| | |l988[ |, it may be too drastic to work with a 
fixed metric for horizontal, low-frequency motion like that associated with the g-modes. This is worth keeping in mind should one want 
to carry out a quantitatively accurate analysis of the problem. 
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Figure 1. Composition gradients for our chosen equations of state (BSk20 from [GurieIy et al.||2013[ | and the particular stellar model 
we focus on throughout the paper. The left-hand panel shows the fractions of protons, electrons and muons versus the total baryon 
density, n, and compares the approximate description used in this work (solid lines) with the analytical fits to the full numerical solution 
obtained by |Potekhin et al.| | |20l'3| l (dashed lines). The right-hand panel shows the core composition for our stellar model, which has 
mass M = IAMq and radius R = 11.66 km. The vertical lines in the right-hand panel represent, respectively, the position where = 0 
{Rfi = 10.53 km) and the crust/core interface {Rcc = 10.82 km). Note that the right-hand panel is using a stretched radial coordinate 
for which the star’s surface is pushed to infinity. Hence, the crust region (indicated as grey in the figure) has been truncated, but this is 
unimportant as we do not consider the dynamics of that region. 


2.3 Composition gradients 


In order to complete the background model and discuss the internal composition of a given star, we obviously need to provide 
an equation of state. In principle, we can assume that this is given in terms of an energy (density) e = e{n,Xx), where 
n = Tin -b np is the total baryon number density and Xx = Hx/n are the various species fractions. Once we have this relation, 


we can readily work out the pressure from (lOl. There are numerous proposed “realistic” equations of state available in the 


literature, so one might think we would be spoilt for choice. Unfortunately (or perhaps fortunately?) this is not the case. 
The problem is that most tabulated equations of state do not provide all the information we need in order to study neutron 
star oscillations. Basically, the tabulated data often assumes that the matter is in equilibrium, whereas we need to be able 
to track how the matter responds when driven out of equilibrium by the fluid motion. The equilibrium configuration must 
satisfy three conditions; 


beta equilibrium: /in = /ip + Me 


lepton balance: 


Mp — Me ) 


and 


charge neutrality: Up = Ue + rip 


( 20 ) 

( 21 ) 


( 22 ) 


but the perturbations may be such that a moving fluid element violates one or more of these conditions. In fact, it is precisely 
this mismatch that leads to the buoyancy that gives rise to the presence of g-modes. In other words, we need to get the physics 
input from a general energy functional that allows us to model matter out of equilibrium. For practical reasons, it would be 
advantageous if the equation of state was given in a parameterised analytic form rather than as a numerical table. We will 
need to work out various partial derivatives and it would be an extremely tricky exercise to do this in a thermodynamically 
consistent way starting from tabulated data (this should be clear from the discussion in Appendix A). For these reasons, we 


base our analysis on the family of equations of state provided by the Brussels-Montreal collaboration (Fantina et al. 2012 


Goriely et al.|[2013 Fantina et al.|[^13 |. Even though we focus on one particular family of models, the parameter space is 
still overwhelming. Hence, we will provide actual results only for a single stellar model. The star we consider is constructed 
from the BSk20 equation of state (Goriely et al.||2013|. A discussion of the general features of this equation of state model 


Fantina et al. 


(20131 and Potekhin et al. (20131. Our particular stellar model has mass 


and stars built from it can be found in 
M = 1.4M0 and radius R = 11.66 km. This should be a fairly “average” neutron star. 

Let us now consider the interior composition of this model. In principle, this problem must be solved numerically. However, 
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in the spirit of the present study it would be preferable to have an analytic approximatiorj^ Fortunately, such an approximation 
is readily obtained if we note that the symmetry energy S{n) (which represents the energy cost in replacing protons with 
neutrons in symmetric nuclear matter) satisfies (see |Chamel||2008 1 

~ he (Syr^nXe)^^^ ~ 45'(n)(l — 2xp) , (23) 

where h is the reduced Planck’s constant and c is the speed of light (which is taken to be unity throughout most of the 
discussion of dynamics, but reinstated in parts of the discussion of the equation of state input), and we assume that the 
electrons are sufficiently relativistic that we can neglect their rest mass. This is a good approximation for a neutron star core 
(and most of the crust). For a typical neutron star, the outer layers of the core only contains neutrons, protons and electrons. 
Then we can use the condition for charge neutrality to replace Xp with Xe on the right-hand side of the equation. Thus we 
obtain 

qT ~ 1 


1 

Srr^n 


45 

he 


1 + 


(24) 


Once the muons appear, the problem becomes a little bit more involved, but in reality (241 remains quite accurate (basically 


since the fraction of charged particles is small). This means that we can solve (211 for the muon fraction. To do this, we need 

1/2 

(25) 


= mp,c 1 -I- 


'cfm j 


T 2 


where kpp, = (ZTi^nXp) 


2/3 


is the Fermi wavenumber for the muons. Thus we arrive at 

1 3/2 


Xu 




m 


— ) (3 


2 l2/3 

TT n 


Y'" - 1 


(26) 


Finally, the proton fraction follows from Xp — Xe + Xp,- As an illustration, the resulting fractions for the BSk20 model are 
shown in the left panel of Figure These results agree well with the parameterised expressions given by [Potekhin et al.| 
(20131, which are based on the numerical solution to the problem. Our approximate fractions are accurate enough for our 


purposes. The right-hand panel of Figure shows the composition of our model star, indicating the crust-core transition 
{Rcc = 10.82 km here) as well as the radius at which the muons first appear {Rp = 10.53 km in our case). 


2.4 Thermal gradients 

The buoyancy of a local fluid element depends not only on the matter composition, it can also be affected by thermal gradients. 
In order to account for this aspect, which should be important during the early stages of a neutron star’s life, we carry out a 
detailed cooling simulation for our model star. 

The temperature evolution of an isolated neutron star is determined by the relativistic equations of energy balance and 
heat flux (see Yakovlev &: Pethick||2004 and Page et al.||2006 for reviews). We solve this system using the method described 
in Ho et al. ( |2012| ). The physics inputs for these equations are the heat capacity, emissivity of neutrinos (which is the primary 
cooling process at ages < 10® yr), and thermal conductivity. The heat capacity is the sum of contributions from the constituent 


particles, also given in Ho et al. (20121. Various neutrino emission processes contribute to the total emissivity, with the primary 


processes being the modified Urea mechanism and Cooper pair formation and breaking (see Ho et aL]|2012 and references 
therein); the BSk20 EoS does not produce neutron stars that undergo neutrino emission via the direct Urea mechanism. 


Thermal conductivities in the core are calculated using results from references given in Ho et al. (20121, while those in the 
crust are calculated using CONDUCTl^ Note that, due to high thermal conductivity the crust and the core are roughly 
isothermal after ~ 10 — 100 yr (see Fig. |^, although there are still slight variations at later times due to the superfluidity. 
Two other important factors affecting neutron star cooling behaviour are the composition of the outer layers (envelope) of the 
neutron star crust and whether the stellar interior is superfluid and/or superconducting. For the former, the envelope serves as 
a heat blanket, with envelopes composed of light elements conducting heat more efficiently than heavy element envelopes. For 
simplicity, we only consider neutron stars with an iron envelope. For the second factor, superfluidity and superconductivity 
suppress heat capacities and some neutrino emission mechanisms, as well as enhancing neutrino emission through the Cooper 
pairing process. 

The initial temperature for our evolutions is taken to be a constant Te* = 10^® K. Figure]^ shows the time evolution of 
temperature at different ages of our model star. Also shown are the critical temperatures for the onset of neutron superfluidity 
in the singlet state in the crust and triplet state in the core and proton superconductivity in the core. Specifically, we use the 


^ Note that the analytic fits for the particle fractions from [Potekhin et al.||2bl3| are not directly useful for the present problem as they 
assume matter at equilibrium and we need to quantify the response to deviations from this equilibrium state. 

® http://www.ioffe.ru/astro/conduct/ 
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-log(l -r/R) 


Figure 2. A sequence of cooling results for particular choices of the superfiuid pairing gaps. The initial temperature for our evolution 
is taken to be a constant Te^ = 10^^ K. We show the critical temperatures for the onset of neutron superfluidity in the singlet state in 
the crust (ns; red solid line) and triplet state in the core (nt; red solid line) and proton superconductivity in the core (p; blue dashed 
line). Specifically, the neutron singlet model is that from [Ainsworth et al.| <|1989||, t he neutron triplet model is taken from [Amundsen &[ 
[0stgaai^ l|1985||, and we use the proton superconductivity model from |Chen et al.| ( |1993| l. The parametrisation of the gaps is the same 
as in [Ho et al. |l |2015| |. Note that the critical temperature for neutron superfluidity in the core is relatively high compared to many other 
models in the literature. We have considered more moderate cases, but focus on this case because it was the only model that led to the 
presence of a brief era of unstable g-modes (see Figure [plater). Note also that thermal gradients are more pronounced in the crust. The 
results in this figure provide input for our discussion of the thermal effects for particular ages of the star {ti in the figure). 


neutron singlet model from Ainsworth et al. ( 1989[ ), the neutron triplet model is taken from Amundsen & 0stgaard (1985 I, and 
we use the proton superconductivity model from Chen et al. ]( p^ . Critical temperatures are calculated using parametrised 
gap models following the prescription from Ho et al. (20151. Note that, as we do not yet have self-consistent calculations of 
the relevant pairing gaps for neutrons and protons for the BSk equations of state (or, indeed, any other proposed equation of 
state), we are “forced” to use this phenomenological prescription. Note that, in the following we do not distinguish between 
singlet and triplet pairing for the neutron superfiuid. In essence, we take the critical temperature to be Ten = max(rcns, Tent)- 
We also assume that T^p > Ten throughout the core (as in Figure]^. 

Our results show that, at very early times, the core cools more rapidly than the crust via stronger neutrino emission, so 
that the crust is generally at higher temperatures. A cooling wave travels from the core to the surface, eventually bringing the 
neutron star to a relaxed, isothermal state. We also see the effect of superfluidity, i.e., faster cooling after neutrons become 
superfiuid in the crust at early times and in the core at later times. These effects are strongest in regions near the critical 
temperatures, in accordance with the expectations from previous work (Gusakov et al.|2004(. 


2.5 Accounting for thermal effects 


Given the results of the recent analysis by Gusakov & Kantor (20131 and the suggestion that thermal effects may lead to 


convective regions in a cooling neutron star, we need to be able to consider both superfluidity and thermal effects. A key 
feature of this problem is that the leptons (electrons and muons) are the main carriers of entropy in a superfiuid neutron star 
core. As long as the local temperature is far below the critical temperature for the onset of superfluidity, it is safe to ignore 
the presence of thermal excitations in the superfiuid. However, these excitations play an important role near the transition. 
In essence, there may always be a local region (near the critical density for the given temperature) where excitations must 


be accounted for. As discussed by Andersson et al. (20131, this effect can be quantified in terms of entrainment between the 


superfiuid and the normal component (entropy). In practice, this leads to the superfiuid component locking to the normal 
component at the transition. A key question concerns to what extent the superfiuid transition is local. If it is, then one may 
be able to avoid a detailed analysis and simply impose suitable junction conditions. If, on the other hand, the interface is 
extended over a macroscopic region, then any study of neutron star dynamics will have to account for the detailed transition. 
As we will account for the presence of thermal excitations, our analysis should shed some light on this issue. 

Our aim is to combine the results of the cooling simulations with a local analysis of the fluid dynamics. When it comes to 
the inclusion of thermal effects, the discussion involves approximations and it is useful to explain what they are. First of all. 
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the cooling simulation is carried out for a stellar model with fixed size, internal composition and so on, by solving the heat 
equation on top of a passive stellar configuration. This approach is standard ( Yakovlev fc Pethick|2004 Page et al.|2006 1 and 
should be adequate as long as one can neglect the effects of thermal pressure compared to the cold degeneracy pressure. As our 
cooling simulation starts at a temperature of 10^° K ~1 MeV and the typical chemical potentials are of order 100 MeV, this 
approximation should be safe. The errors involved are certainly smaller than our level of ignorance regarding the supranuclear 
physics. 

When it comes to quantifying the thermal effects, it is useful to consider the following argument: The total pressure has 
to satisfy 


P = 


x=n,p,e,// 


+ sT' = - I nx/Tx + sT j <!>' = - (p + e) <E>' , 


x=n,p,e,/^ 


(27) 


where we assume that the star is composed of neutrons, protons, electrons and muons (as before) and the entropy density is 


given by s. In a region where the neutrons are superfluid, they must satisfy (191. If we also impose the equilibrium conditions 


(balance between electrons and muons and charge neutrality) we are left with a relation that can be written; 


np^[(Mp+Pe)e*] 


d 


dT^ 

dr 


(28) 


Combining this with (191 we arrive at 


(29) 


dr N - np dr ' 

This relation (also given by Gusakov fc Andersson|2006 ) tells us that a model that satisfies the various chemical equilibrium 
conditions must also be in thermal equilibrium. That is, we must have T°° = constant. This is, of course, not going to be 
true during the early cooling phase, see Figure]^ In order to make the model consistent we would need to add the thermal 
contributions to the relevant chemical potentials. This would then lead to a gradual evolution of the star’s composition during 


the cooling; essentially the tail-end of the de-leptonisation that dominates the hrst 100 s or so of the star’s life (Burrows & 


Lattimer||1986 1. However, as this effect is very small during the phase that we are considering it should be safe to neglect it. 

In quantifying the fluid dynamics we do, however, want to account for thermal effects. We do this by considering the 
equation of state as a zero-temperature model with thermal contributions added perturbatively. Assuming that only the 

where 


and the relevant Fermi temperature is given by 


TT fesTrix 
2 Tfx 


fcsTFx = he (Stt^Wx) 


^ e,p 


1/3 


(30) 


(31) 


(note that, for simplicity, we are ignoring the muon mass here). “Integrating” the definition of the temperature (as the entropy 
chemical potential); 


"-'I 


we see that the thermal energy is 


£th = -^ksT 


UeT n^T 

Tpe Tpu 


and the thermal pressure, which follows from the integrated first law (101 

i. rj, f rieT n„T 
Pth = ( — -h 


(32) 


(33) 


(34) 


Tpe Tp ^_ 

This contribution tends to be (very) small compared to the cold degeneracy pressure, so we only include it whenever it plays 
a leading role in the analysis. 


In the case of a superfluid neutron star core, we will work with S = s/rig (cf. Equation (97l). For this variable, the above 
relations imply that we have 


S' = S 


^ 2 / Tie -b (ne/n^)^^^n(, 

T Tie 3 I Ue -I- (nen2)l/3 


In the low-density region of the core where the muons are yet to appear, this reduces to 

\T' 


S' = S 


T 


n-e 

3ne 


(35) 


(36) 
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3 COMPOSITION G-MODES OF A COLD NORMAL-FLUID STAR 


In order to establish our approach to the problem, and provide useful results for comparison, let us begin our analysis by 
revisiting the standard g-mode problem ( Reisenegger fc Goldreich|1992 I, ignoring superfluidity and thermal effects. That is, 
we consider a stratified star where all fluid constituents move together, in such a way that there is only one fluid four velocity 
to worry about. In this case the equations that govern deviation from a static equilibrium are (i) the perturbed Euler equation 

(e + p) dtSuk = -e'^dkSp - e* {Se + Sp) 9^$ , (37) 


and (ii) the baryon number conservation law 

e~^dtSn + dk + n 


+ <0^ -|— I 5u^ + cot 6 5u^ 
r 


= 0 . 


(38) 


The system of equations is closed once we specify an equation of state. Looking ahead to the problem of main interest, that of 
a superfluid neutron star core, we assume that the star is composed of a mixture of neutrons, protons, electrons and muons. 
Our aim is to establish how variations in composition leads to buoyancy and the emergence of a set of g-modes. 


3.1 The dispersion relation 

We assume that the dynamics is fast enough that the composition of each fluid element can be considered frozerj^ This is 
equivalent to stating that the Lagrangian variation (A) of the particle fractions Xx vanisl|^ This means that we have 

Aa:x = 0 , —> Sx^ = , (39) 

where S denotes Eulerian perturbations and is the Lagrangian displacement vector. In the case we are considering (a static 
background star) we have 

= Cut = e-'^dtC ■ (40) 


where £u represents the Lie derivative along the background four velocity. 

In order to study the dynamics of the star, we assume a harmonic time-dependence and use the standard (polar) 
decomposition of the displacement vector in terms of spherical harmonics We then have (suppressing sums over I 

and m, since the different multipoles decouple for a spherical background star) 






V 


r\ -^rm iat 

OeYi e , 


e = 


V 


r\ -irm icrt 

OcpYi e , 


(41) 

(42) 

(43) 


sin 02 

where W and V are functions of r. Scalar quantities, like the perturbed pressure, are also expanded in spherical harmonics, 
which means that we have Sp = etcetera. In the following, hats are frequently used to denote amplitudes of perturbed 

quantities. These variables are functions of r only. 

In order to make progress on the problem, we need to decide which variables to work with. In this first example, we will 
choose to work with the perturbed pressure, Sp, and the radial component of the displacement vector, W. This means that 
we need to obtain relations for the perturbed number density, 5n, and the perturbed energy density, 5e, from the equation of 
state. As we are assuming that the composition of each fluid element is frozen on the timescale of the motion, it is natural to 
consider the equation of state as a functior0e = e{p, x-Pj. Moreover, it is convenient to work with the enthalpy, w = p -\- e. 
This means that we have 


5w = 


1 / 






(44) 


^ The veracity of this assumption can always be checked a posteriori by comparing the inferred dynamical timescales to the relevant 
dissipative timescales, etcetera. 

^ For later reference it is worth noting how this condition is derived. The conclusion follows from the general result for the Lagrangian 
perturbation of each number density 

Anx = -^ [59^6 + 2V(,C,)] . 

Whatever the fluid motion is, the metric variation Sg^h is the same for all components. If two components are locked, in the sense that 
their displacements are the same, it is easy to show that the Lagrangian perturbation of the ratio of the two number densities must 
vanish. 

® Note that this is not the form we assumed when we discussed the equilibrium configuration. We are using the pressure, p, rather than 
baryon number density, n, as one of the primary variables. This is not necessarily the most practical choice, but it facilitates a direct 
comparison with the work of < 


Gusakov & Kantor 


(2013 I and 


Kantor Gusakov 


120141 
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Here, and throughout the rest of this Section, the sum over x involves p, e and If we define the speed of sound as 

'dp'' 


2 

c, = 


de 


and use (391 we see that (44 \ leads to 


where we have defined (for later convenience) 




1 _ ^ f de 
Tx e V 


y 7^ X. 


P,Xy 


Finally introducing 


we arrive at 




u)= ( 1+^ ]p+-N!w . 


With these various definitions, we also have the radial component of the Euler equation; 

wa We ' ' = p + gw , 

and the angular component; 


2t, -2* 

wa Ve = P■ 

The latter immediately allows us to eliminate V from the analysis. Finally, the baryon number conservation law leads to 


Here we see that we need 


IF'+( — + ^+A')W = 


h Ip 




V - - . 


E p6 ' 


n Fj p 

where 

1 _ p f dn' 
r? n \dpy _ 

(the index b indicates that the adiabatic index is associated with the baryon number n) and 

J_ - i f ^ 

FJ n V dx^ 


y 7^ X. 


Defining 


we have 


A? = 




where At = X^x^?- Introducing also 


n P A 

- = —=^ - WAb 

n pFl 


r2 _ P + 1) 2«y^i) 
*-!) —-5—e 11 


we have the two equations 


and 




pT 


P + » + ^) P = y - N^] W 


where 


2 - £2(<S>-A) ,.2 


K' = -e 

w 


Nt . 


(45) 

(46) 

(47) 

(48) 

(49) 

(50) 

(51) 

0 

(52) 

(53) 

(54) 

(55) 

(56) 

(57) 

(58) 

(59) 

(60) 

(61) 


In order to analyse the local dynamics of the problem we now make the plane-wave assumption, i.e. assume that all 
perturbations depend on position as e**’’ where kr 3> 1, as the focus is on the local dynamics. We also note that we can 
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introduce integrating factors on the left-hand side of each equatiorj^ However, it is straightforward to demonstrate (see 
Andersson fc Comer|200T| that these integrating factors do not affect the final dispersion relation. The upshot of this is that 
we can effectively “ignore” the second term of each equation. Thus, we arrive at the two relations (ignoring the overall 
factor, as this should not cause any confusion) 


ikW 


p 2 A-<r 
- r p 

r\a^ p 


(64) 


Combining these we have 


where 


ikp = e — A/)?) ^ . 


p V 




fc = e* ^k . 


(65) 


( 66 ) 

(67) 


This should be the complete dispersion relation, without approximations other than the plane-wave assumption. We see 
that there are two sets of solutions; representing travelling waves in two limits. As we need k^ > 0, we must have either 
> Afc (generally expecting <C jC^) or cr^ < A/)? < Cl- The first set of solutions correspond to high-frequency 
sound waves given by; 


2 


a 




w 


p l{l + 1) 2*p6 _ p2 _ Id + 2*^ / dp \ 

w r'^ ^ ^ ’’ r'^ ^ w \ dn) ^ 


( 68 ) 


The second approximation is valid whenever kr I, an assumption we will use to simplify the multi-fluid problem later. 
Meanwhile, for low frequencies we have (composition) g-modes with frequency; 


1 - 


Kyl+l) 




.r2 £ 2(*-A) 

: A4 = -g—e 
w 


E 


Tx 


(69) 


This result shows how the gradients of the different particle fractions provide the buoyancy that leads to the g-modes. We 
also see that the modes will be unstable if the leading gradient increases outwards in the star (provided Tx > 0). This is the 
familiar criterion for the onset of convection. Of course, considering the actual composition gradients from Figure we see 
why convection is not expected to take place in cold neutron stars. The composition stratification is stable (although it is 
worth noting that the there is a region where the electron contribution has a destabilising effect in this particular example). 

Figure]^ shows numerical results for our model star (for 1 = 2). We see that the muons make the dominant contribution 
to the g-modes throughout the part of the core where they are present. The upshot is that the g-mode frequency may be up 
to an order of magnitude higher than in pure npe matter | |Kantor k, Gusakov||2014| |. We also see that (for this rather low 
value of 1) the sound waves may, in fact, have lower frequency than the g-modes in the outer part of the core. This basically 
shows that, if we want to assume that 3> A/)?, as in the above discussion, then we need to choose a larger value of 1. It is 
useful to keep this in mind later. 


3.2 Thermodynamical relations and results 

Most previous g-mode calculations have been carried out using a simplified approach to the equation of state, where the 
stratihcation is accounted for by making the adiabatic index for the perturbations different from that for the background 
(which has often been taken to be a simple polytrope). This approach is fine as long as one is mainly interested in qualitative 
properties. It is also, although perhaps to a lesser extent, natural given the difficulty to extract the required information for 
tabulated realistic equations of state. As we are aiming for a higher level of realism, we need to address the full problem. This 
means that, in order to obtain the results shown in Figure]^ we had to work out various partial derivatives for our equation 
of state. As this is an important step, and one has to work things out in a thermodynamically consistent way, let us consider 
this problem in more detail. For clarity, we focus on the region where all four particle species are present. 


For completeness; the integrating factors we need are, for Equation ||59[l; 


nr exp 


/ 


A — / Abdr 


(62) 


-fg/c 




and for Equation ( |60| l; 


exp 


(63) 
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-log(l -r/R) 


Figure 3. Propagation diagram for a cold, normal fluid, neutron star core. We show the frequencies u = a'/27r for local sound waves 
(represented by for I = 2 obtained from l |68| l) and the g-modes (given by obtained from | |69| l). The results show that the muons 
make the dominant contribution to the g-modes throughout the part of the core where they are present. 


First of all, we have the variation of the energy density (taken to be in the original form we assumed when we determined 
the background configuration) 

p + £ 


de = = 


-dn + n{fip — fin)dxp + npedxe + nppdxp 


x=n,p,e,;z 


where we have used the (integrated) first law (lOl. From this relation, it is easy to see that 


2 

c, = 


dp 

dn 


Flence, we see that our sound-wave solutions can be written 


2 l{l + 1) 24 > 

o- = A = „ e 


[ 1 ^ 5 ] 


_ P-rb 
— 1 
w 




(70) 


(71) 


(72) 


which means that our notation makes sense. It is also worth noting that the factor in the square brackets can be evaluated 
without reference to a specific stellar model. 

This comparison was quite straightforward. However, in order to evaluate, for example, ( |47[ ) we need to hold the pressure 
p fixed in the variation. This means that we need to ensure that dp = 0. We can recast this constraint into a relation that 
gives dn in terms of the variations in the particle fractions; 

-1 




dn 


v=n 


dp 


dx„ 


dx^ 


y 7^ X . 


(73) 


Using this in (7^ we can calculate the different partial derivatives required to determine the individual Fx. These can then 
be used in (|^ I to obtain the (local) g-mode frequencies for a given stellar model. Adding the different terms together we get 
(note: the variables that are held fixed are not stated explicitly here) 


E w = 


1 1 
F 

in i f> 


_1 _]_ 

,F„, F, 


p + e 


dp 

&n 


dp 

dxp 


dp 

dxg 


dp 

dxu 


dp 

dxg 


(74) 


In the first equality we have assumed that the background is charge neutral (22l. In the second equality we have imposed the 
background beta-equilibrium (201 as well as the balance between muons and electrons (211. Of course, as we have already 
discussed, the background composition can be expressed in terms of the baryon number density, so in fact we have Xp = Xp{n), 

leans that the fi 

dp dp \ dx 


and similar for all other background quantities, which means that the final expression we need in (691 is 

-1 r 


E x'x w I dp 

P.. n \ i9r) 


dn 


dxp dXf 


dn 


dp 

dXn 


dp 

dxe 


dXp 

dn 


(75) 


x=p,e,/i 

In order to quantify the g-mode frequencies, and obtain results like those in Figure we obviously need to build a stellar 
model. However, it is worth noting that, if the focus is on the relative contributions of the proton and muon gradients, then 
all we need to do is compare the two terms in the square bracket in (751. This can be done without reference to a particular 
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star as it only requires the partial derivatives and the known composition as a function of n. Specifically, we first note that 

-1 / I \ -1 


p = -gw 


(76) 


which means that we have 


2 2 2(‘i»-A) 

(7 = g e ^ 



(77) 


To quantify the terms in the square bracket we only need the pressure and composition of matter at equilibrium. We will 
not provide an explicit example of this comparison, because the information it would contain is already clear from Figure]^ 
However, it is still useful to keep this possible comparison in mind because it would provide an immediate and simple way of 
checking the importance of different contributions to the buoyancy for a variety of equations of state. 


4 THE G-MODES OF A SUPERFLUID NEUTRON STAR CORE 

Moving on to the issue that motivated this investigation, we consider how the problem changes when the neutrons become 
superfluid and we also account for finite temperature effects. As we have already discussed, the neutral superfluid component 
may drift relative to the charged components, which remain locked. This adds a degree of freedom to the problem. It also 
means that we need to adjust our notion of composition gradients. Since the neutrons now form a separate component, that 
can move out of the way if it is pushed, they no longer contribute to the buoyancy that the other components experience. The 
upshot of this is that, if we ignore thermal effects, then there is no local support for g-modes in regions where only neutrons, 
protons and electrons are present. The neutrons do not play a role because they are superfluid and charge neutrality requires 
the electron and proton fractions to be equal. Hence, there is no buoyancy and no g-modes ( Lee||1995 Andersson & Comer 
20011. Interestingly, this means that entropy gradients will be key to determining the nature of the g-modes in this region. 


This was pointed out by Gusakov & Kantor (20131, who also argued that these gradients may be such that the resulting 


modes are unstable. Of course, at the density where muons appear the situation changes again. Beyond this density there will 
again be a composition gradient as the ratio of muons to protons and electrons will vary ( [Kantor fc Gusakov||2014| ). This is 
likely to overwhelm the thermal effects and should serve to suppress any instability. Our aim in the following is to test the 
veracity of these expectations. 

In a region where the neutrons are superfluid, the locked components still satisfy conservations laws of the form (381 
[replacing n with rip, rie, or s = Us]. Meanwhile, we have to account for the perturbed drift velocity, in the conservation 
law for the neutrons. This means that we have 


= 0 , 


(78) 


e~*9t(5nn -I- dk ^rin (Su'^ + -1- >1?' -|—j {5u" + Sv^) -I- cot 0 

Adding this to the equation for the protons we have the conservation law for baryons; 

e~'^dt5n + dk -\- nn(5u*^j -I- ^A' -I- -I- {nSu^ + UnSv’’) + cot 6 (nSu^ + UnSv^'^ 

If we ignore entrainment (we will account for this effect later) then the (perturbed) momentum equation for the superfluid 
neutrons takes the very simple form; 


= 0 . 


(79) 


dtSpl = dkSpt ■ 


That is, we have 


Pn {dtSuk + dtSvk) = -e^dk&Un - Spr^e^dk^ = -dk , 


(80) 


(81) 


where the last identity holds since we are working in the Cowling approximation where the metric is held fixed. The final 
equation we need represents total momentum conservation. As the neutrons are allowed to drift relative to the other component 
this relation is different from ( |37| |. We now have 

wdtSuk + UnUndtSvk = —e^dkSp — dk^Sw , (82) 

where we have (again) used the enthalphy w = p + e. 

In this problem we need to keep track of two Lagrangian displacements. Following previous work on the problem, we take 
these to be and , defined such that (cf. the velocity dependence in (781 and @); 

-I- XnSv °‘, (83) 

and 


— $ ^ a c a I c a 

e dtp = ou -1-01! . 


( 84 ) 
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This means that (81|-(82| become, respectively 




k d TJk 2^ O r 2'i' C 

— c^n = —e Okop — ge dw ^ 


dt^ 


where we have defined 


d Pk 2 ' 1 'q c * 

~~Qp2~ ~~ ^ 1 


5pl = 


5p,n 

/^n 


We have also used the background equation (191. Finally, we have defined 

Oc = - « W{T = 0) + 


sT 


and 


(85) 

( 86 ) 

(87) 


( 88 ) 


Xj,W - UnAtn Xn . ■, Xr 

— — = - yW - WMn) ~ — 

Xp Xp Xp 


sT 


(89) 


The last equalities assume the cold equilibrium conditions for the chemical potentials, which is somewhat inconsistent as 
we wonld then have to have T = Q. However, the purpose here is to argue that <C Qc which means that we can often 
neglect this contribution. Having said that, we will retain this quantity throughout the calculation because this will facilitate 
a comparison with the case where the entrainment is accounted for later. 


4.1 The dispersion relation (no entrainment) 


In order to derive the dispersion relation for the superfluid problem, we (again) assume time dependence and decompose 
the variables in spherical harmonics. The displacement is decomposed as in (41l-(43l, while ? 7 “ takes the same form but 
with Wn and Vn replacing W and V. As in the analysis of the normal-fluid problem in Section 3, we use hats to denote the 
amplitudes of various perturbed scalar quantities. 

The two scalar conservation laws then lead to 


W' + 


W^n + 


n' 2 
— -b - -b A' 
n r 

^ + - + A' 


_nn r 

while the radial components of the momentum equations become; 


w=^-^Pv-^, 

n 

Wn = - — 


rin 


and the angular parts give; 


Combining these we have 


cr^ (acW — OnlFn) e ^^=drP + gw, 
= drPl , 


cr^ (q-cH - OnK) e 


2-,t - 2 €> 

cr 146 


f> + OnAn 2« 

V = -^ 6 


(90) 

(91) 

(92) 

(93) 

(94) 

(95) 

(96) 


and we see that we can remove both V and 14 from the analysis at this point. 

Let us now assume that we focus on the part of the star’s core where muons are present, and that we choose to work with 
p, pu and S = s/ue and M = Ufijne as our primary variable^ These are the same variables as Gusakov & Kantor (20131 


used, which facilitates a direct comparison. We assume that the perturbations are adiabatic, which means that 


A5 = 0 


6S = -CsS' = 


nnW4 — nW 


S', 


(97) 


where is the radial component of the displacement associated with (and can be obtained by inverting (831 and ( |84[ )). 
If the muons and electrons are also locked together (as we assume), then the relation for SM is obtained simply by replacing 
S with M in the above result. In fact, as these two quantities enter the equations that follow in a similar way it is convenient 
to introduce a label X which can be either S ox M (and also use Y 7 b X). 


® Note that S is defined as the entropy per electron. This makes sense because the electrons and muons carry the entropy when the 
neutrons are superfiuid and the protons are superconducting. We also need the ratio to be defined between two co-moving components 
in order to show that the Lagrangian variation of this quantity vanishes, cf. Section 3.1. 
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Making use of thermodynamics, we have for the enthalpy: 


5w = 5e + 5p=[l + ^]5p+ - -NI - W) , 


(98) 

y 

where the term that ultimately leads to the buoyancy takes the form 

Here, and in the following, we use bars over the adiabatic indices (eg. Fx) to make a distinction from the analogous quantities 
in the normal-fluid case discussed earlier. It is important to make this clear as we use different variables and different quantities 
are held fixed in the partial derivatives. 

As in the case of a system where all components are locked together, the buoyancy depends on a balance of terms. In the 
present case we have thermal stratification and the muon gradient. It is worth noting that overall charge neutrality implies 
that 


M' = 


- 1 


( 100 ) 


from which it is easy to see that the second contribution to vanishes in the region where the muons are absent (as it must). 

( 101 ) 


In (981 we have used 


-2 

Cc - 




(as before, the bar is used to make a distinction between the superfluid and the normal-fluid problem) 

1 _ Un f de 




p,X 


(where the index /r on F refers to the chemical potential rather than the muons) and 

Jx - i 

Tx ^ e ' 

Similarly, the neutron number density perturbation is now given by 

SUn _ 1 , 1 r * , W-nlFn - uW 

- ^ tvT-^ -’ 

n-n Ff p F“ rip 


with 


/ flOn 
ff Tin V dp 


( 102 ) 


(103) 


(104) 


(105) 


Pii,X 


(where the n index indicates that these adiabatic indices are associated with the variation of the neutron number density rin 


1 


F" rin V dp 


dn^ 


p,X 


and, for later convenience, we have defined 




F" 

X 


with 


1 _ J_ ^ dUn ^ 

Finally, we have the baryon number perturbation; 


rx \ dX J 


P,l^n,Y 


Sn 1 Sp 1 c * rinVFn — nW . 

— — -f ^—Ab , 

n ri p F* rip 


where we have used 


1 

= ^ / 

^dn\ 

F[ 

n ( 

ydp)^ 

1 

/^n 

/ dn 

pF 

n 



(106) 

(107) 

(108) 

(109) 

( 110 ) 
( 111 ) 


p,X 


and 


. X' 

X 


( 112 ) 
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with 


1 1 f dn 


X 


\dX 


P,IJ.n,Y 


As in the normal-fluid case, it is natural to introduce; 


r'2 _ P + 1) 2Spi) 
■c-i) —-5—e i 1 


(113) 


(114) 


(Uc r'= 

Note that, all quantities labelled C (with various indices) contain the “centrifugal” l{l -\- 1) factor. They do not, in the general 
setting, directly correspond to solutions to the final dispersion relation. That was particular to the normal-fluid problem. 


Let us now manipulate equation (90l; 

W' + 


+ ? + A' _ 


r2 2 - 

w = - + 


f\a2 p 


l{l -b 1) 

ac<7^ 


Yb 


f-^n -A-h 


Similarly, equation (911 becomes 


VT' + 


— + - + A' + — A 




An CT 


Ip An 
+ - W 

Ffp Xp 


where we have defined 


2 _ Id + 1) 2«pn 

•'-n “ ® ^ P ■ 


Next consider the radial part of ( |92[ ); 

drP -l-p(l-l-4)p = ace 


-2(*-A) 


(cr^ — A/’c ) W + i XuHc — CF 


2 P^r. 


Wn 


- df^Pn , 
^ P 


where we have defined (cf. the analogous definition in the cold, normal fluid problem, Eq. fsTf) 


etc 


Nf: 


Finally, the momentum equation for the neutrons is simple: 

drpu = , 


(115) 

(116) 

(117) 

(118) 

(119) 

( 120 ) 


The four equations (1151, (1161, (1181 and (1201 provide us with all the information we need to determine the dispersion 
relation in this more complex setting. As in the single-fluid problem, we can introduce integrating factors to simplify the 
equation^ and yet again the upshot of this is that we can effectively ignore the terms in brackets on the left-hand sides 
of ( [TT^ , ( |116[ ) and ( |118[ |. Making this simplification and focussing on plane-wave solutions (i.e. taking the perturbations to 
depend on the radial coordinate as but not changing the notation for the amplitudes, as there should be no risk of 

confusion) we have 


r2 2 - 

Ak — O’ T) 

ikp = 


Id 1) OnC^ 

J.2 

J^P 
r? P 


+ 


CTcCr 

A 


pi, 

P2 


U-nAb 

Mn- Wn , 


W, 


Xp 


(P — Afp) W + i XnJVp — a 


Wp 


- fXPn, 
^ P 


ik/l*=PWne . 


(122) 

(123) 

(124) 

(125) 


The steps required to obtain the dispersion relation are straightforward. First we remove W and Wn from the equations 
to get two coupled equations for p and /!„. Finally, combining those equations, we arrive at the dispersion relation 


p2 2 

, 2 . ( 2 ^/-2^ ‘-b “ “ 

p ^ ‘ F^a^ 




pn 
^ M 


jXnk ik / 2 \r2\ 

H—— AA ) 


pry 


/' 2\ ^ I 4 2(-i>—A) 

/i 


ik^ 


,r2 2 . gS 

acXnXc - ancr + ^ 

k r^j 


ik 


An cl- > 


pFJ XpP rjer^ 


(126) 


The integrating factors for equations ( |115[ |, i prel i and ( |ll8| l are, respectively, given by 


n r exp 


A- 

f^dr] 


J Xc 


Tin exp 

A- 

f —Andr 

, and exp 

[/»(.+i)*i 



J Pic 


Li V clJ \ 


( 121 ) 
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Figure 4. The two sound-wave solutions that follow from the simplified dispersion relation | |128[ l for the superfluid problem (and 1 = 2)^ 
assuming that the star is sufficiently cold that the entire core is superfiuid. Note that the solutions show an avoided crossing in the outer 
core of the star. We also show (as a dashed line) the speed of sound from the normal-fluid problem (L^ from | |68[ |). The results show 
that one of the two sound speeds in the superfiuid problem tends to be close to the “normal sound”, but a switch-over takes place at the 
avoided crossing. 


where we have introduced 


r2 fS& + 1) 2$ 

2 ^ 


(127) 


and used k = ke^~^, as before. It is worth noting that the dispersion relation is now a cubic in a^, which means that we 
expect to find three (more or less distinct) classes of waves. In addition to the acoustic modes and the g-modes from the 
single-fluid problem we should have a set of superfiuid modes ([Andersson fc Comer|200T|. 


4.2 Results for a cold superfiuid core 


The result (1261 provides the complete dispersion relation, without approximations other than the plane-wave assumption. Let 


us now add to this the assumption that we are dealing with short-wavelength/large multipole waves, such that 1 <C fcr <C 1 (as 


before). After comparing the different part of (1261 and keeping only the dominant terms, we arrive at the simplified relation 
for high-frequency waves 




where 


X = 


r b 
1 ^ 

■rn pb 

1 /i 

The sound waves follow from this quadratic in cr^. Formally, we get 


2 

Cl,2 — 


(1 ± V>)£g + (1 T V>)^g 

2(1-X) 


where 


= 


1-f 


4x^6^g 


1/2 


(128) 


(129) 


(130) 


(131) 


This result shows that the sound waves are typically some combination of and Only in the particular case when x = 0 


do these quantities themselves represent solutions to the problem. This agrees with the results of Andersson & Comer (20011. 


In practice, we also find that one has to be a little bit careful when evaluating these expression. For our chosen equation 
of state, the quantity FJ) changes sign in the outer core and this could cause numerical difficulties. The actual sound-wave 
solutions are, however, regular. Results for our model star (and I = 2) are shown in Figure The frequencies generally 
decrease outwards in the star, as expected, and we also see an avoided crossing; a common feature of this kind of problem. 
Comparing to the normal fluid case, we find that one of the two solutions in the superfiuid problem tends to be close to the 
“usually assumed” sound speed. However, the assignation changes at the avoid crossing, which is an interesting observation. 
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Figure 5. A comparison of the two contributions to the superfiuid g-modes. The results show that the muons contribution always 
dominates over the thermal part even if we artificially inflate the temperature to 10^^ K, an order of magnitude above our initial 
value for the cooling simulations. 


Turning to the g-modes, the low-frequency solutions to (1261 are approximately given by 

ge ^ 2 («>-a) 


= M^ = — 

dcXp 


S' M' 


— (J s + (T M 


(132) 


where the coefficients Ts and Tm are discussed in Appendix A. In this expression, the muon contribution to the buoyancy 
can be written (making use of Oc = w) 
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while the entropy contribution takes the form 
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where we have used 


Ts 


(133) 


(134) 


(135) 


and the coefficient A is obtained by combining (301 with the results from Appendix A2. 

Making use of these results, we can carry out two comparisons. First, we can establish that the muons contribution aj^,i 
always dominates over the thermal part (j|. This is shown in Figure where we have artificially inflated the temperature 
to 10^^ K, an order of magnitude above our initial value for the cooling simulations (in this comparison the temperature is 
assumed to be uniform, so we only include the first term in (1341). Of course, the thermal contribution has no competition in 
the part of the region where the muons are not present, so we still need to consider it. The second comparison is between the 
superfluid g-modes and the normal fluid results from Section 3. This comparison is made in Figure]^ assuming that the star 
is cold. This means that there is no support for g-modes until the muons appear. In the region where the muons are present, 
the g-mode frequency is found to be significantly higher in the superfluid case. This is simply understood from the fact that 
the neutron have been decoupled, which means that the mass involved in the wave differs by a factor of (roughly) Xp. That 
this is, indeed, the difference between the two cases can be shown by taking a closer look at the two expressions for Afc- 

Let us now combine our results with the cooling data from Section 2.4. That is, let us trace how the g-modes evolve as 
the star cools and the core becomes superfluid. If we ignore the entrainment (including thermal excitations) then the results 
of this exercise should be easy to understand. In a region where the fluid is above the superfluid transition temperature, the 
normal-fluid results from Section 3 should apply, and in superfluid regions we have to use the results from this Section. How 
this works out in practice is shown in FigureThe results clearly show how the superfluid region grows as the star cools, in 
accordance with the results from Figure 

Finally, let us consider the suggestion from Gusakov & Kantor (20131 that the superfluid g-modes may become unstable 
and trigger a convective phase in a young neutron star. From the results we have discussed so far it is clear that, if such an 
instability is to occur it has to be located in the outer part of the crust where the muons are absent. In essence, we need the 
overall sign of the terms in the bracket of (1341 to be negative for this instability to be present. Considering this issue for our 
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Figure 6. A comparison of the superfluid g-modes (from | |132| l) and the normal fluid results (from l |69[ l), assuming that the star is cold 
so that thermal effects can be ignored. In the superfluid case there is no support for g-modes until the muons appear. In the region where 
the muons are present, the g-mode frequency is found to be significantly higher in the superfluid case. This is simply understood from 
the fact that the neutrons have been decoupled, which means that the mass involved in the wave differs by a factor of (roughly) Xp. 


chosen equation of state and the single stellar model we are focussing on, we have carefully checked the various temperature 
distributions extracted from our cooling simulation. From this data we have found one single case where an instability is 
present. This case is illustrated in Figure]^ We see that the instability is weak ((j| only becomes marginally negative) and 
strongly localised, meaning that only very short wavelength motion would actually be unstable. The instability is present 
after about a day of thermal evolution (which is long enough that the artificial isothermal initial temperature distribution 
has filtered through the system, as this happens on a timescale of seconds), but it is gone well before the system is one year 
old. This seems to suggest that this convective instability is unlikely to play much of a role in the evolution of a neutron 
star. Having said that, it is conceptually interesting and one should keep in mind that we have only considered one particular 
model. 


4.3 Entrainment 


If we want to account for the relevant physics of a neutron star core, it is imperative that we include the entrainment effect in 
our model. Entrainment is known to play an important role in superfluid dynamics | Mendell||1991[ Prix et al.|2002 Andersson 
et al.|[Ml2 |. This is natural as it encodes how easy (or hard) it is for a superfluid component to move relative to the other 
parts of the system. It is a non-dissipative effect, which essentially represents the effective mass of the superfluid component. 
There are two, formally different but de facto equivalent, ways of incorporating the entrainment. The first approach encodes 
the effect in a mass-density matrix pxy which designates how much one component tends to flow along with the other(s) 
I Mendeli]||1991 1. An alternative, which is more natural from the relativistic point-of-view ( Andersson fc Comer|[2007 1, is to 
account for the effect through the momentum of each component. In this approach, the entrainment leads to the momentum 
of a given component pj no longer being proportional to the corresponding flow nj, but rather a linear combination of all the 
flows that are entrained. 

Under the conditions that prevail in a cold neutron star core, the main cause of entrainment is the strong interaction. 
As a superfiuid neutron moves, a virtual cloud of entrained protons tries to move with it (and vice versa). The effect does 
not depend (at least not much) on the temperature. In the following we quantify this strong interaction entrainment by a 
coefficient In a finite-temperature superfluid, one also has to account for the presence of thermal excitations. As argued 
by Andersson et al. (20131, the dynamical role of the excitations can be accounted for in terms of entrainment between 


the superfluid and the heat ( Lopez-Monsalvo fc Andersson||2011 Andersson fc Lopez-Monsalvo||2011 1. In the following, we 
represent this effect by a coefficient A"'*. The upshot is that the momentum of the superfluid neutrons in our mixture is given 
by (see |Andersson et al.||2013| for a detailed discussion) 

/x“ = B’^nl + A'^'^nl + A’^'^Sa , (136) 

where the three coefficients S", and follow from the equation of state. It is worth noting that we have assumed that 
the protons, electrons and muons all flow together. Projecting the neutron momentum along the neutron flow, we see that 
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Figure 7. A sequence of snapshots corresponding to temperature distributions extracted at three particular times from Figure]^ The 
results show how the superfluid region grows as the star cools, and how we have to replace the normal-fluid g-mode frequency from 
Section 3 with the superfluid results from Section 4. The speciflc cooling data relates to Iq = 50 yr, = 60 yr and fg = 70 yr, 
respectively. The corresponding temperature profiles are shown in Figure]^ 


the corresponding chemical potential is given by 


Mn = + sA^ 


(137) 


This expression is the key to adding the entrainment to our previous analysis of the g-mode problem. This extension turns 
out to be surprisingly straightforward. After consulting the discussion from the previous Section, we simply introduce a new 
parameter 


P = = 1 - T + sA"^) . 

f-^n 


(138) 


As discussed in Appendix B, there are two key limits to consider. First of all, in the limit of zero entrainment we have /3 = 1. 
Meanwhile, /3 will diverge as T —^ Ten, as we approach the critical density for onset of neutron superfluidity (see Appendix B, 


especially Eq. (B13l). 


The analysis leading to the dispersion relation now proceeds along the same steps as in the case without entrainment. 
We only have to account for the fact that the momenta are linear combinations of the two fluxes and (carefully) redefine some 
of the quantities involved. We will also need slightly different combinations of the various thermodynamical derivatives. The 
logic of the analysis is, however, the same as before. Hence, we relegate the details to Appendix B and focus on the results 
here. The main question is how the entrainment impacts on the local wave propagation. 


Let us first consider the sound waves. In the high-frequency limit, we find that the dispersion relation (B35l reduces to 

{Cl - P) {Cl - ^ ^ 0 , (139) 

where 


and 
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P - Xn 
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In these expressions, Cl takes the same form as in (11141), but we now have 
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P - Xn 


Meanwhile, 
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(143) 


In general, we have the two anticipated sets of sound waves. However, the situation changes as we approach the critical 
temperature for the onset of superfluidity. In this limit, we have 

(144) 


lim ac = w , 

T-fT^a 


and 
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Figure 8. An illustration of the presence of a region of unstable g-modes. This example corresponds to the single one of our extracted 
temperature distributions for which an instability is present, at ta = 10“^ yr (cf. Figure]^. We see that the instability is weak and 
strongly localised, meaning that only very short wavelength motion would actually be unstable. 


This means that 
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(147) 


We learn that one set of sound waves disappears as we cross the superfluid transition. This is, of course, as expected. The 
“second sound” is not supported in normal matter. As far as we are aware, this is the first time that this limit has been 
demonstrated explicitly. In particular, our analysis highlights the key role that the thermal entrainment plays in the problem. 


A closer analysis of the remaining sound wave, represented by (1461, shows that it is close to, but not identical to, the normal 
sound (f 1 is not identical to TJ). The typical difference is at the few percent level, depending on where the superfluid transition 
takes place in the core. In other words, a small discontinuity remains in the local sound-wave frequency as we cross the critical 
density between superfluid and normal regions. 


Turning to the low-frequency regime, we find that (B351 leads to the presence of g-modes with frequency 

2 ^ 
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where 
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which leads to, as we approach the superfluid transition; 
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as T —^ Tci 


(148) 


(149) 


(150) 


As in the case of the sound waves, this result does not limit to the normal-fluid g-mode frequency at the superfluid transition. 
In fact, the discontinuity in this case is significant. The superfluid results starts out significantly below the normal-fluid 
g-mode frequency at the superfluid transition temperature (which accords with the results of Kantor fc Gusakov||2014 1 and 
then rises towards the cold superfluid result as the star ages. Making use of our cooling data, we obtain the results shown in 
Figure]^ This figure illustrates the effect that thermal excitations (accounted for through the entrainment between superfluid 
and entropy) has on the transition between different layers in the star. Note that we have not included the entrainment due 
to the strong interaction as is it essentially temperature independent (and it is a small effect). This demonstration shows that 
it is important to account for the thermal effects when considering the dynamics of superfluid neutron stars, and emphasises 


the need for more detailed modelling building on, for example, Gusakov & Andersson (20061; Kantor & Gusakov (20111; 


Chugunov & Gusakov ( 2011|; Andersson et al. (20131; Gualtieri et al. (20141. 
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Figure 9. An illustration of the effect that thermal excitations (accounted for through the entrainment between superfluid and entropy) 
has on the transition between different layers in the star. The superfiuid result starts out significantly below the normal-fluid g-mode 
frequency at the superfiuid transition temperature (which accords with the results of |Kantor &c Gusakov|2014[ | and then rises towards 
the cold superfiuid result as the star ages. The results relate to stars aged tg = 10^ yr, tio = 10^ yr, tn = 10^ yr and ti 2 = 10® yr, 
respectively, cf. the cooling curves in Figure]^ A key take-home message is that the oscillation spectrum may still be evolving even for 
rather mature stars. 


5 CONCLUDING REMARKS 


In this paper we have investigated the local dynamics of a realistic neutron star core, accounting for composition gradients, 
superfluidity and thermal effects. Our main focus was on the gravity g-modes, which are supported by composition stratification 
and thermal gradients, but we also provided results relevant for the sound waves of the system. We have derived the detailed 
equations that govern the problem, paying particular attention to how the physics that is encoded in the equation of state 
enters and distinguishing between normal and superfiuid regions. The analysis highlighted a number of issues that need to be 
kept in mind whenever equation of state data is compiled from nuclear physics calculations for use in this kind of analysis. 
We provided explicit results for a particular stellar model and a specihc equation of state (BSk20 from Goriely et al.|[20T3 l, 
making use of accurate cooling simulations to show how the local wave spectrum evolves as the star ages from adolescence to 
maturity. Our results conhrm the expectation that the composition gradient is dominated by the muons beyond the density at 
which they hrst appear ( Kantor fc Gusakov|2014 1. At lower densities, where there are no muons, the support for the g-modes 
is entirely thermal once the star cools below the superfiuid transition. We confirm the recent suggestion that the g-modes in 
this region may be unstable ( Gusakov fc Kantor||2013 1, but our results indicate that this instability will be weak and would 
only be present for a short period of the star’s life (less than a year). A novel technical aspect of our analysis is associated 
with the thermal excitations at hnite temperatures which are accounted for in terms of entrainment between the entropy and 
the superfiuid component. Once the thermal entrainment is accounted for we see that the oscillation spectrum may still be 
evolving even for rather mature stars. Finally, we discussed the difference between the normal sound waves and the second 
sound that is present only in the superfiuid regions. 

This work extends the state-of-the-art in several ways. Most importantly, we have shown how detailed equations of state 
data can be implemented to study neutron star dynamics at a serious level of realism. This poses two main challenges for 
the future. First of all, we need to build on this work to construct global models for dynamical neutron stars, as required for 
detailed asteroseismology studies. Secondly, we need a continued dialogue with nuclear physics experts both to expand the set 
of equation of state models that provide all the information required for this kind of study, and also to complete the models 
we have used in this work (e.g. concerning consistent superfiuid pairing gaps). 
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APPENDIX A: REQUIRED THERMODYNAMICAL RELATIONS 

Let us now take a closer look at how we evaluate the different thermodynamical partial derivatives we need. This issue is 
not trivial as our perturbative analysis made use of variables that require a somewhat unusual description of the equation of 
state. Basically, it is important to keep track of what is being held fixed when the various adiabatic indices are evaluated. As 
a first step, which takes us a bit closer to where we need to be, let us assume that the equation of state is given in terms of 
the variables n = + rip, Xe = rie/n, M = n^lrie. and S = s/n-e- As our main interest is in the g-modes, let us focus on the 

quantities we need to determine them. The other thermodynamical quantities can be obtained in an analogous fashion. 

In terms of the chosen variables, the energy variation can be written 

de =+ nSTdxe + nXeTdS , (Al) 

n 

(where we have imposed the relevant equilibrium conditions to simplify the prefactors). 


Al The muon contribution 


Let us, first of all, work out the contribution due to the muon gradient. This means assuming that S is held fixed, so we take 
dS = 0 in the following. We want to work out 

As this involves holding the pressure constant, we use the definition of the pressure to get 

dn = —aidxe — a 2 dM , (A3) 

where 
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(A4) 


(here and in the following we need to keep in mind that the other variables are to be held fixed when the partial derivatives 
are evaluated) and 
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The last equality follows since n and Xe are held fixed, which means that dM = {l/xe)dXfj,. Next, the requirement that /in is 
constant leads to the relation 


dxe = —a^dM 


where 
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By combining these results, we arrive at the expression we need 




(A6) 


(A7) 


(AS) 


A2 The entropy contribution 

Turning to the thermal contribution, we instead hold M fixed and work out 

’ de \ 


For obvious reasons, the results are similar to those for the muons. Holding the pressure fixed leads to 

dn = —h\dXe — b2dS , 

where 

h ^{dpyWdp^ 
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(the only difference from the previous case is that M is held fixed rather than S) and 
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Next, keeping /in fixed leads to the relation 
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By combining these results, we arrive at 
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In the model we consider in this paper, the thermal contributions are added perturbatively to a cold background star. 
Moreover, it tends to be the case that the thermal effects can be neglected compared to the cold parameters. For example, 
for the range of temperatures we consider, the thermal pressure can always be neglected compared to the cold degeneracy 
pressure throughout the star’s core. For our model, we find that the thermal pressure leads to 
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This provides us with the information required to evaluate the 62 coefficient above. Meanwhile, for the other partial derivatives, 
we can neglect the thermal pressure. Writing the total pressure as p = po + Pth we then have 


dp\ ^ (dpo\ 


dn 


dn ) 


and 


Turning to the neutron chemical potential, we first of all have 
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The other partial derivatives we need can be accurately obtained from the cold equation of state. 
Combining the results, we have the final relation 
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where we note that the overall factor of T implies that the thermal support for the g-modes will decrease as the star cools. 
This is, obviously, as expected. 

In practice, we can use 
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where the terms in the brackets can be obtained from the cold equation of state. 


A3 Other thermodynamical relations 


In addition to the thermodynamical derivatives required to determine the g-modes, we need a more complete set if we are to 
solve the full dispersion relation. These involve; 
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It is worth keeping in mind that, for the model we are using, we have 
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APPENDIX B: ACCOUNTING FOR ENTRAINMENT 


B1 Entrainment coefficients 


Following Andersson et al. 120131, we encode the entrainment between the superfluid neutrons and thermal excitations in terms 
of a coefficient which can be obtained from the results of Gusakov & Haensel (20051. They provide a (non-relativistic) 
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mass density matrix, given by; 


Pnn = (1 - /n) Pnn , (Bl) 

Ppp = (1 “ /p) Ppp 5 (B2) 

Pnp = (1 — /n) (1 — /p) Pnp , (B3) 

where /x are temperature-dependent functions that tend to zero in the T —>■ 0 limit and approach unity when T —>■ Tex- The 

particular form for these functions is given by Gnedin & Yakovlev (19951. In the above expression, the Pxy quantities depend 
only weakly on the temperature. Hence, we take them to be given by the zero temperature expressions of the mass density 
matrix in the following. Expressing the results in terms of the effective proton mass m* we then have (Andersson et al.|2013( 
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where px = mn^ is the non-relativistic mass density and the approximations are valid when rip/rin <C 1. 


The translation to our relativistic model has already been discussed by Andersson et al. (20131. It proceeds in two steps. 


First we identify the correspondence with the leading order Newtonian terms; 
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where we have defined 


(l-/n)7^ 

Ti- = Pnnppp — (1 — /n) (1 — fp) Pnp ■ 


(B7) 

(B8) 
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In the relativistic case, we also need to add the Newtonian chemical potential p,^ to rinH". This term is (within our approxi¬ 
mations) temperature independent. We then have 

Pn = -I- Upk'^^ + = m + . 


The last equality shows that the neutron chemical potential is not affected by the presence of thermal excitations. 
In the cold limit, we can use (B4l-(B6l (not the approximations, though!) to show that 




Tl-p 771 Pnp 


(Bll) 


(B12) 


Pn detp 

as expected. Meanwhile, as the critical temperature is approached we have the leading order behaviour 

m PnPpp ^ 

Pn (1 /n)^ 

As we will see, the dynamics of the system is controlled by the divergence of this quantity as T —>■ Ten. 

Finally, it is worth pointing out that, if we ignore the thermal effects then the strong interaction contributes to the 
entrainment exactly as in a cold superfluid neutron star. We would have 


A:“p = — fm - ; 


(B14) 


Although important, this effect is relatively small. As our main focus is on the thermal contributions we have not accounted 
for it in the results shown in Figure 


B2 Dispersion relation (with entrainment) 


When we turn to the perturbation problem, we see that the entrainment requires a slight redefinition of the displacement 
vectors. Specifically, we have 
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The radial components of the momentum equations are then given by 

cr^ [a^W — anVKn) = drP + gw , 

= drp*^ , 

while the angular parts are 

cr^ {aP/ - a„Vn) = p , 

2-,r -2* --* 

cr VnC = /in • 

In the previous expressions, we have defined the following quantities: 
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XxiW 7ln/^n 

/3-Xn " ■ 

The two conservation laws lead to 
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Tin 


where we have used equation (B23I to introduce n/n in equation!B24 1 . The expansions of Sw, Srin/n^ and 5n/n are similar 

;o Xx, as 

IT^ x' . 


to the zero entrainment case, the only difference is in the terms related to x^, as we now have: 
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Explicitly, we have: 
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With these definitions equation (B241 can be re-written as 

P - Xn [ P - I 
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where we have defined 
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(B32) 

(B33) 

(B34) 


The steps in the derivation remain exactly as in the case without entrainment. Carrying out the required algebra, we 
arrive at the final dispersion relation 
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where we have used the definitions; 


and 
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AC = ^Ar2e2(*-A) ^ 
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Results obtained from this dispersion relation are discussed in the main body of the paper. The key point to note from our 
analysis is that, even though /3 diverges as we approach the critical temperature, all equations (and hence the dynamics!) 
remain regular. 
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